function image_compression()
% IMAGE_COMPRESSION 基于SVD的图像压缩
% 
% 演示SVD在图像压缩中的应用

fprintf('=== 基于SVD的图像压缩 ===\n\n');

%% 1. 创建测试图像
fprintf('1. 创建测试图像\n');

% 创建多种类型的测试图像
img_size = 128;

% 图像1：几何图案
[X, Y] = meshgrid(1:img_size, 1:img_size);
img1 = sin(2*pi*X/32) .* cos(2*pi*Y/32) + 0.5*sin(4*pi*X/64) + 0.3*cos(6*pi*Y/48);
img1 = (img1 - min(img1(:))) / (max(img1(:)) - min(img1(:))) * 255;

% 图像2：随机纹理
rng(42);
img2 = conv2(randn(img_size), ones(5)/25, 'same') * 50 + 128;
img2 = max(0, min(255, img2));

% 图像3：分块图像
img3 = zeros(img_size);
block_size = 16;
for i = 1:block_size:img_size
    for j = 1:block_size:img_size
        img3(i:min(i+block_size-1,img_size), j:min(j+block_size-1,img_size)) = ...
            randi([0, 255]);
    end
end

images = {img1, img2, img3};
image_names = {'几何图案', '随机纹理', '分块图像'};

fprintf('创建了%d个测试图像，尺寸: %d×%d\n', length(images), img_size, img_size);

%% 2. SVD图像压缩分析
fprintf('\n2. SVD图像压缩分析\n');

for img_idx = 1:length(images)
    img = images{img_idx};
    img_name = image_names{img_idx};
    
    fprintf('\n--- %s ---\n', img_name);
    
    % SVD分解
    [U, S, V] = svd(img);
    singular_values = diag(S);
    
    fprintf('奇异值统计:\n');
    fprintf('最大奇异值: %.2f\n', singular_values(1));
    fprintf('最小奇异值: %.2e\n', singular_values(end));
    fprintf('条件数: %.2e\n', singular_values(1) / singular_values(end));
    
    % 奇异值分布分析
    total_energy = sum(singular_values.^2);
    cumulative_energy = cumsum(singular_values.^2) / total_energy;
    
    fprintf('能量分布:\n');
    energy_thresholds = [0.9, 0.95, 0.99, 0.999];
    for threshold = energy_thresholds
        rank_needed = find(cumulative_energy >= threshold, 1);
        compression_ratio = rank_needed / img_size;
        fprintf('保留%.1f%%能量需要%d个奇异值 (压缩比: %.3f)\n', ...
                100*threshold, rank_needed, compression_ratio);
    end
    
    % 不同压缩级别的效果
    ranks_to_test = [1, 5, 10, 20, 40, 80];
    fprintf('\n压缩效果分析:\n');
    fprintf('秩    压缩比    PSNR(dB)    SSIM      MSE\n');
    fprintf('--    ------    --------    ----      ---\n');
    
    for r = ranks_to_test
        if r <= min(size(img))
            % 压缩重构
            img_compressed = U(:, 1:r) * S(1:r, 1:r) * V(:, 1:r)';
            
            % 压缩比
            original_size = numel(img);
            compressed_size = r * (size(img, 1) + size(img, 2) + 1);
            compression_ratio = compressed_size / original_size;
            
            % 图像质量评估
            mse = mean((img(:) - img_compressed(:)).^2);
            psnr = 10 * log10(255^2 / mse);
            ssim_val = compute_ssim(img, img_compressed);
            
            fprintf('%2d    %6.3f    %8.2f    %6.4f    %7.2f\n', ...
                    r, compression_ratio, psnr, ssim_val, mse);
        end
    end
end

%% 3. 彩色图像压缩
fprintf('\n3. 彩色图像压缩\n');

% 创建简单的彩色图像
rgb_size = 64;
[X_rgb, Y_rgb] = meshgrid(1:rgb_size, 1:rgb_size);

% RGB三通道
img_r = sin(2*pi*X_rgb/16) .* cos(2*pi*Y_rgb/16) * 127 + 128;
img_g = cos(2*pi*X_rgb/20) .* sin(2*pi*Y_rgb/20) * 127 + 128;
img_b = sin(2*pi*X_rgb/12) + cos(2*pi*Y_rgb/12) * 127 + 128;

img_rgb = cat(3, img_r, img_g, img_b);
img_rgb = max(0, min(255, img_rgb));

fprintf('彩色图像尺寸: %d×%d×%d\n', size(img_rgb));

% 分别对每个通道进行SVD压缩
target_rank = 10;
img_rgb_compressed = zeros(size(img_rgb));

channel_names = {'红', '绿', '蓝'};
fprintf('\n各通道压缩分析:\n');
fprintf('通道    原始能量    压缩能量    能量保持率\n');
fprintf('----    --------    --------    ----------\n');

for c = 1:3
    channel = img_rgb(:, :, c);
    [U_c, S_c, V_c] = svd(channel);
    
    % 压缩
    img_rgb_compressed(:, :, c) = U_c(:, 1:target_rank) * ...
                                  S_c(1:target_rank, 1:target_rank) * ...
                                  V_c(:, 1:target_rank)';
    
    % 能量分析
    original_energy = sum(diag(S_c).^2);
    compressed_energy = sum(diag(S_c(1:target_rank, 1:target_rank)).^2);
    energy_ratio = compressed_energy / original_energy;
    
    fprintf('%s    %8.0f    %8.0f    %10.3f\n', ...
            channel_names{c}, original_energy, compressed_energy, energy_ratio);
end

% 整体压缩效果
total_compression_ratio = (target_rank * (rgb_size + rgb_size + 1) * 3) / numel(img_rgb);
fprintf('\n彩色图像压缩结果:\n');
fprintf('压缩比: %.3f\n', total_compression_ratio);

% 计算彩色图像的PSNR
mse_rgb = mean((img_rgb(:) - img_rgb_compressed(:)).^2);
psnr_rgb = 10 * log10(255^2 / mse_rgb);
fprintf('PSNR: %.2f dB\n', psnr_rgb);

%% 4. 渐进式图像传输
fprintf('\n4. 渐进式图像传输\n');

% 选择一个测试图像
test_img = images{1};
[U_prog, S_prog, V_prog] = svd(test_img);

fprintf('模拟渐进式传输过程:\n');

% 渐进式传输的秩序列
progressive_ranks = [1, 2, 5, 10, 20, 40, 80];
fprintf('传输阶段    累积秩    数据量(KB)    PSNR(dB)    视觉质量\n');
fprintf('--------    ------    ----------    --------    --------\n');

cumulative_data = 0;
for stage = 1:length(progressive_ranks)
    current_rank = progressive_ranks(stage);
    
    % 累积数据量（简化计算）
    if stage == 1
        data_increment = current_rank * (size(test_img, 1) + size(test_img, 2) + 1);
    else
        prev_rank = progressive_ranks(stage-1);
        data_increment = (current_rank - prev_rank) * (size(test_img, 1) + size(test_img, 2) + 1);
    end
    cumulative_data = cumulative_data + data_increment;
    data_kb = cumulative_data * 8 / 1024;  % 假设双精度
    
    % 重构图像
    img_progressive = U_prog(:, 1:current_rank) * S_prog(1:current_rank, 1:current_rank) * ...
                      V_prog(:, 1:current_rank)';
    
    % 质量评估
    mse_prog = mean((test_img(:) - img_progressive(:)).^2);
    psnr_prog = 10 * log10(255^2 / mse_prog);
    
    % 视觉质量评估（简化）
    if psnr_prog < 20
        visual_quality = '差';
    elseif psnr_prog < 30
        visual_quality = '一般';
    elseif psnr_prog < 40
        visual_quality = '良好';
    else
        visual_quality = '优秀';
    end
    
    fprintf('%8d    %6d    %10.1f    %8.2f    %s\n', ...
            stage, current_rank, data_kb, psnr_prog, visual_quality);
end

%% 5. 自适应压缩
fprintf('\n5. 自适应压缩\n');

% 根据图像内容自适应选择压缩参数
fprintf('自适应压缩策略:\n');

for img_idx = 1:length(images)
    img = images{img_idx};
    img_name = image_names{img_idx};
    
    [U_adapt, S_adapt, V_adapt] = svd(img);
    singular_values_adapt = diag(S_adapt);
    
    % 分析图像复杂度
    % 方法1：基于奇异值衰减速度
    sv_ratios = singular_values_adapt(1:end-1) ./ singular_values_adapt(2:end);
    avg_decay_rate = mean(sv_ratios);
    
    % 方法2：基于有效秩
    sv_normalized = singular_values_adapt / singular_values_adapt(1);
    effective_rank = sum(sv_normalized > 0.01);
    
    % 方法3：基于能量集中度
    energy_90 = find(cumsum(singular_values_adapt.^2) / sum(singular_values_adapt.^2) >= 0.9, 1);
    
    fprintf('\n%s复杂度分析:\n', img_name);
    fprintf('平均衰减率: %.2f\n', avg_decay_rate);
    fprintf('有效秩 (>1%%): %d\n', effective_rank);
    fprintf('90%%能量秩: %d\n', energy_90);
    
    % 自适应选择压缩参数
    if avg_decay_rate > 5  % 快速衰减，简单图像
        recommended_rank = min(20, energy_90);
        complexity = '简单';
    elseif avg_decay_rate > 2  % 中等衰减
        recommended_rank = min(40, effective_rank);
        complexity = '中等';
    else  % 慢衰减，复杂图像
        recommended_rank = min(60, effective_rank);
        complexity = '复杂';
    end
    
    fprintf('图像复杂度: %s\n', complexity);
    fprintf('推荐压缩秩: %d\n', recommended_rank);
    
    % 验证推荐参数
    img_auto_compressed = U_adapt(:, 1:recommended_rank) * ...
                          S_adapt(1:recommended_rank, 1:recommended_rank) * ...
                          V_adapt(:, 1:recommended_rank)';
    
    auto_compression_ratio = (recommended_rank * (size(img, 1) + size(img, 2) + 1)) / numel(img);
    auto_mse = mean((img(:) - img_auto_compressed(:)).^2);
    auto_psnr = 10 * log10(255^2 / auto_mse);
    
    fprintf('自适应压缩结果: 压缩比=%.3f, PSNR=%.2f dB\n', ...
            auto_compression_ratio, auto_psnr);
end

%% 6. 压缩性能比较
fprintf('\n6. 压缩性能比较\n');

% 与其他压缩方法的概念性比较
test_img_perf = images{1};
[U_perf, S_perf, V_perf] = svd(test_img_perf);

fprintf('不同压缩方法的理论比较:\n');
fprintf('方法        压缩比    质量    计算复杂度    适用场景\n');
fprintf('--------    ------    ----    ----------    --------\n');
fprintf('SVD低秩     中等      高      高            科学图像\n');
fprintf('JPEG        高        中      低            自然图像\n');
fprintf('PNG         低        无损    中            计算机图形\n');
fprintf('小波变换    高        高      中            多分辨率\n');

% SVD压缩的优缺点分析
fprintf('\nSVD图像压缩特点:\n');
fprintf('优点:\n');
fprintf('• 数学上最优的低秩近似\n');
fprintf('• 适合具有低秩结构的图像\n');
fprintf('• 支持渐进式传输\n');
fprintf('• 压缩比可精确控制\n');
fprintf('• 理论基础扎实\n\n');

fprintf('缺点:\n');
fprintf('• 计算复杂度高 O(mn²)\n');
fprintf('• 不适合高频细节丰富的图像\n');
fprintf('• 压缩比通常不如专用算法\n');
fprintf('• 需要存储浮点数\n');
fprintf('• 块效应可能出现\n\n');

%% 7. 实际应用考虑
fprintf('\n7. 实际应用考虑\n');

% 计算复杂度分析
img_sizes = [64, 128, 256, 512];
fprintf('不同图像尺寸的计算复杂度:\n');
fprintf('尺寸      SVD时间(相对)    存储需求(相对)\n');
fprintf('----      -------------    --------------\n');

for size_val = img_sizes
    % 复杂度估算（相对于64×64）
    base_size = 64;
    time_complexity = (size_val/base_size)^3;  % O(n³)
    storage_complexity = (size_val/base_size)^2;  % O(n²)
    
    fprintf('%4d      %13.1f    %14.1f\n', ...
            size_val, time_complexity, storage_complexity);
end

% 实际应用建议
fprintf('\n实际应用建议:\n');
fprintf('1. 预处理:\n');
fprintf('   • 图像分块处理大图像\n');
fprintf('   • 预先分析图像复杂度\n');
fprintf('   • 考虑颜色空间转换\n\n');

fprintf('2. 压缩策略:\n');
fprintf('   • 根据应用需求选择压缩比\n');
fprintf('   • 结合感知质量评估\n');
fprintf('   • 考虑渐进式传输需求\n\n');

fprintf('3. 优化技术:\n');
fprintf('   • 使用快速SVD算法\n');
fprintf('   • 并行计算加速\n');
fprintf('   • 量化减少存储需求\n\n');

fprintf('4. 适用场景:\n');
fprintf('   • 科学可视化图像\n');
fprintf('   • 医学影像压缩\n');
fprintf('   • 卫星遥感图像\n');
fprintf('   • 教育演示用途\n\n');

%% 总结
fprintf('=== SVD图像压缩总结 ===\n');
fprintf('SVD提供了图像压缩的数学最优解，特别适合:\n');
fprintf('• 具有低秩结构的图像\n');
fprintf('• 需要精确控制压缩比的场合\n');
fprintf('• 渐进式图像传输\n');
fprintf('• 科学和工程应用\n\n');

fprintf('在实际应用中，需要权衡:\n');
fprintf('• 压缩效果 vs 计算成本\n');
fprintf('• 图像质量 vs 存储需求\n');
fprintf('• 通用性 vs 专用性\n\n');

end

% 辅助函数
function ssim_val = compute_ssim(img1, img2)
% 简化的SSIM计算
mu1 = mean(img1(:));
mu2 = mean(img2(:));
sigma1_sq = var(img1(:));
sigma2_sq = var(img2(:));
sigma12 = mean((img1(:) - mu1) .* (img2(:) - mu2));

c1 = (0.01 * 255)^2;
c2 = (0.03 * 255)^2;

ssim_val = ((2*mu1*mu2 + c1) * (2*sigma12 + c2)) / ...
           ((mu1^2 + mu2^2 + c1) * (sigma1_sq + sigma2_sq + c2));
end